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ABSTRACT 

We re-analyzed the HerschelfPACS data of a sample of 55 brown dwarfs (BDs) and very low mass stars with spectral types ranging 
from M5.5 to LO. We investigated the dependence of disk structure on the mass of the central object in the substellar regime based 
on a homogeneous analysis of Herschel data from flux density measurements to spectral energy distribution (SED) modeling. The 
broadband SEDs were compiled by adding previous photometry at shorter wavelengths and (sub-)millimeter data. We performed 
detailed SED analysis for the 46 targets that show infrared (IR) excess emission using radiative transfer models and evaluated the 
constraints on the disk parameters through Bayesian inference. A systematic comparison between the derived disk properties and 
those of sun-like stars shows that the disk flaring of BDs and very low mass stars is generally smaller than that of their higher mass 
counterparts, the disk mass is orders of magnitude lower than the typical value found in T Tauri stars, and the disk scale heights 
are comparable in both sun-like stars and BDs. We further divided our sample into an early-type brown dwarf (ETBD) group and a 
late-type brown dwarf (LTBD) group by using spectral type (= M8) as the border criterion. We systematically compared the modeling 
results from Bayesian analysis between these two groups, and found the trends of flaring index as a function of spectral type also 
present in the substellar regime. The spectral type independence of the scale height is also seen between high-mass and very low-mass 
BDs. However, both the ETBD and LTBD groups feature a similar median disk mass of 1 x 10“^ Mq and no clear trend is visible in 
the distribution, probably due to the uncertainty in translating the far-IR photometry into disk mass, the detection bias and the age 
difference among the sample. Unlike previous studies, our analysis is completely homogeneous in Herschel/PACS data reduction and 
modeling with a statistically significant sample. Therefore, we present evidence of stellar-mass-dependent disk structure down to the 
substellar mass regime, which is important for planet formation models. 
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1. Introduction 


How brown dwarfs (BDs) form is one of the main open 
questions in the field of star formation and remains a sub¬ 
ject of debate although several scenarios have been proposed, 
for instance, a scaled down version of star formation pro¬ 
cesses, gravitational instabilities in disks and ejection of the 
stellar embryo (e.g.. iReinurth & Clarke 200 it ^te et al.l l2003t 


IStamatellos & WhitworthI l201 it Chabrier et al.f 20141) . Similar 

to their higher mass counterparts, such as T Tauri stars, young 
BDs are shown to have circumstellar disks, producing sub¬ 
stantial excess emission at wavelengths ranging from infrared 
(IR) to (sub)millimeter (e.g., Liu et al. l2003|^ Scholz_e^i] 2006; 
Ba yo^tVl 2(H2t iHarvev et ^ 2012a . 201 J^cci et al.l 2014 


Liu et al.l ' 20151) . Disks are also found around very faint ob¬ 
jects with masses down to the plan etary regime, for ex - 
ample Cha 110913- 773444 (~8 Mjim , iLuhman et al.l l^05ah . 


LOri 156 (~23 [Bavo et aLll201^. an d OTS44 (~12Mjup, 
iLuhman et al.ll2005bt 1j oergens et al.ll2013h . The phenomena of 


^Jup? 


mass accretion and outflow, wh ich are common by-products of 
the st ar formation process (e.g.. lAudard et al.]l2014l:lLrank et al.l 
1201 4 for reviews), have also been detected in young BDs (e.g. 


* Herschel is an ESA space observatory with science instruments 
provided by European-led Principal Investigator consortia and with im¬ 
portant participation from NASA. 


4ohant^^t_^ 2005 iPhan-Bao et al.ll2()()^ iRigliaco et al.ll201^ 


Joergens et al.l 2012 ) ev en down to the plan etary mass regime 
(e.g., Bavo et alT 2012t Lloergens et al.l l2()T^ . The dust evolu¬ 
tion in BD disks appears to follow a similar manner (i.e., grain 
growth, settling and crystallization) to that in T Tauri disks, al¬ 
though observations suggest different timescales of d ust process¬ 
ing in disks around different stellar mass hos t s (e.g., ^ai^^D 
120051: iBouv et all 120081: iPascucci et al.ll200^ iRiaz et al.ll2012h . 
These observations show that BDs resemble hydrogen-burning 
stars in many aspects during their early evolution, implying that 
they may form through the canonical star formation processes. 


Characterizing the physical properties of disks plays a cru¬ 
cial role in understanding the formation mech anisms of BDs 
fe.g.. IReinurth & ClarkduOOTt iBate et al.ll200^ and planets. In 
addition, thorough comparisons between disk properties like the 
mass, flaring and scale height in different stellar mass regimes 
are very useful and important. So far, many works have been 
done in this direction. Early disk studies in s ubstellar reg i me are 
mainly based on Spitzer data. Lor instance, ISzucs et51 (1201 fll) 
investigated the IRAC and MIPS24yum photometry of ~200 
stars in the Chamaeleon I star formation region and found that 
disks around cooler objects are generally flatter than the case 
of earlier type stars. The Herschel space telescope has unprece¬ 
dented sensitivity and angular resolution in the far-IR, enabling 
the detection of many faint disks at this wavelength domain for 
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the first time. The Herschel data on one hand improve the de¬ 
termination of disk properties, and on the other ha nd provide an 
alternat ive way to estimate the disk mass of BDs. iHarvev et al.l 
(l2012ah observed a large sample of low-mass stars and BDs 
with HerschellPACS and found that the disk masses of low-mass 
stars and BDs extend to well below typical values found in T 
Tauri stars. The lower disk mass around cooler stars as revealed 
by far-IR measurements is consistent with the strong correla¬ 
tion be tween the stellar and disk masses given by [Andrews et al.l 
(I 2 OI 3 I) from a millimeter survey of Class II sources in the 
Taurus molecular cloud. Several recent works obtained similar 
results, although the sample size and assumptions used to ana¬ 
lyze the Herschel data differ (e.g., Alves de^)liyeira_e^.l l201^ 


ISnezzi et'aDl2013t lOlofsson et al 


Ml iLiu et akllMl). Note 


that the above disk comparisons were mostly conducted between 
two stellar mass bins, i.e., the sun-like stars and the BDs includ¬ 
ing very low mass stars. These works have yielded some evi¬ 
dence of spectral-type-dependent disk structure in the low stellar 
mass regime. 

To date, however, there are seldom studies to investigate the 
dependency of disk properties on the spectral type (or the mass) 
of the central object within the substellar regime. It is quite un¬ 
clear that whether the observed trend of disk properties as a func¬ 
tion of spectral type is also valid down to very low-mass BDs and 
even toward planetary mass objects. This kind of study needs 
great care because the disk parameters are not direct observables 
and determining them needs additional assumptions like the disk 
model and dust opacity used in the spectral energy distribution 
(SED) analysis. Moreover, the Herschel data nowadays is better 
understood and this knowledge leads to difference with reported 
fluxe s and uncertainties published at the early stages of the mis¬ 
sion (iBalog et al.ll2014l) . which will in turn induce discrepancy 
in the value of parameters used to describe the disk structure. 
Considering these factors, it is not appropriate to analyze the 
results collected from the literatures that reduce and model the 
Herschel data in diffe rent ways. This may be another reason why 
lOlofsson et al.l (1201 3l) did not hnd a correlation between the stel¬ 
lar and disk masses in the substellar mass regime from a com¬ 
parison between their results and those in the literature. 

In this work, we re-analyzed the HerschellPACS data of a 
sample of 55 BDs and very low mass stars, and ht their broad¬ 
band SEDs using self-consistent radiative transfer models with 
as few free parameters as possible. We estimated the constraints 
on key disk parameters (i.e., flaring index, scale height and disk 
mass) through Bayesian analysis. We extensively discussed the 
trends on these parameters as a function of spectral type. Since 
our analysis is homogeneous from data reduction to modeling 
within a statistically signihcant sample, we present the hrst com¬ 
prehensive study of stellar-mass-dependent disk structure down 
to the planetary mass regime. 

The structure of this paper is as follows: we describe our 
sample in the following section. In Sect. |2 we delineate the ob¬ 
servations and data reduction. A detailed description of the mod¬ 
eling is presented in Sect. |4] We discuss our results in Sect. |5l 
followed by a brief summary in Sect. |6] 


identihcation work and ambiguity of the results. Evidence for 
circumstellar dust emission as revealed by thermal IR measure¬ 
ments is a crucial criterion of our sample selection. 

The sample of this paper was selected based on the study 
by iHarvev et al.l (12012al) on one hand and a search for the low¬ 
est m ass objects in the H erschellPACS archive on the other 
hand. IHarvev et al.l (12012ah conducted the hrst comprehensive 
HerschellPACS program to measure far-IR emission from young 
BDs. We re-perform the analysis for the objects in their survey 
that satisfy the above criteria. Eor each of these targets, we will 
provide the modeling results that were not reported in their work. 
In addition, we include the lowest-mass BDs (with SpT M8 and 
later) observed with PACS from the Herschel archive to probe 
differences in disk properties between high-mass and low-mass 
BDs. Eurthermore, we also add several targets without clear ex¬ 
cess at near- and mid-IR, which are thought to have experienced 
signihcantly inside-out disk evolution, but the cool and outer re¬ 
gions of the disk, if survived, are sensitive to Herschel obser¬ 
vations. We note that including all substellar objects from the 
HerschellPACS archive in our analysis - although desirable - was 
not possible due to limited computational power. The sample 
presented here is therefore biased towards very low-mass sub¬ 
stellar objects. 

Since the SpT is a good proxy for stellar mass in the low 
stellar mass regime, our sample can be generally divided into 
two groups: (1) a group of 30 early-type brown dwarfs (ETBDs) 
with spetral type earlier than M8, (2) a group of 25 late-type 
brown dwarfs (LTBDs) with spetral type M8 and later. With this 
kind of division, we can investigate the dependence of the disk 
structure on the mass of the central object in the substellar mass 
regime, which will yield important clues to understand the for¬ 
mation mechanism of BDs and planets. 

Table [T] summarizes the properties of our sample, includ¬ 
ing target name, coordinates and SpT. Our targets are located 
in eight different star formation regions with ages ranging from 
~l-3Myr (e.g., Taurus and Ophiuchus) to ~10Myr (e.g., TW 
Hydrae association and Upper Scorpius), see Table |2]for a sum¬ 
mary of the cloud ages and distances. The effects of a mixed 
sample age on the results are discussed in detail in Sect. |5] Note 
that for our sample we cannot evaluate the disk frequency of 
BDs at far-IR due to the incompleteness of sample selection and 
the detection bias in various clouds at different distances to the 
Earth. 


2. Sample 

Our sample consists of 55 spectroscopically conhrmed targets 
in nearby star formation regions. The spectral types (SpTs) of 
the sample are in the range from M5.5 to LO. Thus most of 
the sources have subste llar masses accordin g to theoretical evo¬ 
lutionary models (e.g.. iBaraffe et al.]|2003l) . We do not include 
any BD candidate lacking spectroscopic conhrmation to avoid 


2 



































Table 1. Source summary. 


ID 

Name 

RA 

(J2000) 

DEC 

(J2000) 

ObsIDs 

Fy 70 yum 

[mJy] 

Fy 160yum 

[mJy] 

SpT 

Cloud 

Ref 

SD 

Multiplicity 

1 

SSSPM1102 

11 02 09.83 

-34 30 35.5 

1342221849/50 

9.80 ± 1.03 

8.02 ± 1.10 

M8.5 

TWA 

L 2 



2 

2M1207 

12 07 33.47 

-39 32 54.0 

1342202557/58 

7.24 ± 1.03 

< 17 

M8.0 

TWA 

3,2 


binary 

3 

2MJ160859 

16 08 59.54 

-38 56 27.6 

1342227823/24 

7.4 ± 2.0 

<35 

M8.0 

Lupin 

4,2 



4 

GY92 3 (ISO032) 

16 26 21.90 

-24 44 39.8 

1342238816/17 

49.75 ± 5.58 

< 647 

M8.0 

Oph 

5,6 



5 

GY92 264 

16 27 26.58 

-24 25 54.4 

1342238816/17 

33.26 ± 5.45 

< 1489 

M8.0 

Oph 

5,6 



6 

GY92 310 

16 27 38.63 

-24 38 39.2 

1342238816/17 

120.28 ± 5.45 

< 1032 

M8.5 

Oph 

7,6 

Y1 


7 

CFHTWIR96 

16 27 40.84 

-24 29 00.7 

1342227841/42 

< 193 

< 1192 

M8.25 

Oph 

8, 2 



8 

OTS44 

11 10 09.34 

-76 32 17.9 

1342238816/17 

4.63 ± 1.03 

< 211 

M9.5 

Chal 

9, 10 



9 

2MJ111145 

11 11 45.34 

-76 36 50.5 

1342223476/77 

2.84 ± 1.03 

<33 

M8.0 

Chal 

11, 2 



10 

J125758.7-770120 

12 57 58.70 

-77 01 19.5 

1342212708/09 

- 

< 243 

M9.0 

Chall 

12, 13 



11 

KPN0 6 

04 30 07.24 

+26 08 20.8 

1342227011/12 

<4 

<47 

M9.0 

Tau 

14, 2 



12 

KPN0 7 

04 30 57.19 

+25 56 39.5 

1342227999/8000 

2.61 ± 1.03 

<69 

M8.25 

Tau 

14, 2 



13 

J04574903+3015195* 

04 57 49.03 

+30 15 19.5 

1342243478/79 

< 8 

< 14 

M9.25 

Tau 

15, 16 



14 

J04354526+2737130* 

04 35 45.26 

+27 37 13.1 

1342243466/67 

< 11 

< 14 

M9.2 

Tau 

17, 16 



15 

J04334291+2526470* 

04 33 42.92 

+25 26 47.0 

1342243084/85 

< 7 

< 20 

M8.75 

Tau 

17, 16 



16 

J04335245+2612548 

04 33 52.46 

+26 12 54.9 

1342243080/81 

8.00 ± 2.32 

<26 

M8.5 

Tau 

18, 16 

Y1 


17 

J04290068+2755033 

04 29 00.68 

+27 55 03.4 

1342243066/67 

<9 

<36 

M8.25 

Tau 

17, 16 



18 

KPNO 4* 

04 27 28.00 

+26 12 05.3 

1342242064/65 

<9 

< 13 

LO 

Tau 

15, 16 



19 

J04215450+2652315* 

04 21 54.51 

+26 52 31.5 

1342242033/34 

< 7 

<23 

M8.5 

Tau 

17, 16 



20 

J04263055+2443558 

04 26 30.55 

+24 43 55.9 

1342242023/24 

8.00 ± 2.22 

<33 

M8.75 

Tau 

17, 16 

Y1 


21 

J04274538+2357243* 

04 27 45.38 

+23 57 24.3 

1342242021/22 

<9 

<23 

M8.25 

Tau 

17, 16 



22 

J04302365+2359129 

04 30 23.66 

+23 59 13.0 

1342242009/10 

< 7 

< 13 

M8.5 

Tau 

17, 16 



23 

J04311907+2335047* 

04 31 19.07 

+23 35 04.7 

1342242007/08 

< 8 

< 15 

M8.0 

Tau 

19, 16 



24 

KPNO 9 

04 35 51.43 

+22 49 12.0 

1342241979/80 

< 12 

< 16 

M8.5 

Tau 

20, 16 



25 

J04361030+2159364 

04 36 10.31 

+21 59 36.5 

1342241975/76 

< 8 

< 16 

M8.5 

Tau 

17, 16 



26 

J04325119+1730092* 

04 32 51.20 

+ 17 30 09.2 

1342241957/58 

<9 

< 19 

M8.25 

Tau 

21, 16 



27 

KPNO 1* 

04 15 14.71 

+28 00 09.6 

1342241914/15 

<6 

< 19 

M8.75 

Tau 

15, 16 



28 

KPNO 12 

04 19 01.27 

+28 02 48.7 

1342241888/89 

<6 

<49 

M9.25 

Tau 

15, 16 



29 

J04221332+1934392 

04 22 13.32 

+ 19 34 39.2 

1342241878/79 

< 10 

<24 

M8.0 

Tau 

17, 16 


binary 

30 

J04414489+2301513 

04 41 44.90 

+23 01 51.4 

1342240750/51 

< 12 

< 17 

M8.5 

Tau 

17, 16 

Y1 

binary 

31 

ISO 138 

11 08 18.51 

-77 30 40.8 

1342218699/700 

4.52 ± 1.02 

<89 

M6.5 

Chal 

22, 2 



32 

ISO 217 

11 09 52.16 

-76 39 12.8 

1342223474/75 

10.47 ± 1.10 

< 507 

M6.25 

Chal 

23, 2 



33 

CHAHA6 

11 08 39.52 

-77 34 16.7 

1342223486/87 

14.50 ± 1.10 

<54 

M5.75 

Chal 

23, 2 



34 

2MJ110703 

11 07 03.69 

-77 24 30.7 

1342223480/81 

28.48 ± 1.02 

< 447 

M7.5 

Chal 

24, 2 

Yl, Y2 


35 

CHAHA9 

11 07 18.61 

-77 32 51.7 

1342223482/83 

15.30 ± 1.01 

< 220 

M5.5 

Chal 

23, 2 

Y1 


36 

ISO 252 

11 10 41.42 

-77 20 48.1 

1342223478/79 

8.81 ± 1.10 

< 122 

M6.0 

Chal 

23, 2 



37 

ISO 165 

11 08 54.97 

-76 32 41.1 

1342223470/71 

45.46 ± 1.02 

31.89 + 7.03 

M5.5 

Chal 

23, 2 



38 

CHAU 125 8 

12 58 06.76 

-77 09 09.5 

1342224206/07 

5.06 ± 1.01 

< 116 

M6.0 

Chall 

25, 2 



39 

CHAU 1308 

13 08 27.14 

-77 43 23.3 

1342224204/05 

6.43 ± 1.01 

<37 

M6.0 

Chall 

25, 2 



40 

CFHTTAU12 

04 33 09.46 

+22 46 48.7 

1342227013/14 

< 7 

<54 

M6.5 

Tau 

17, 2 
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Table 1. continued. 




ID 

Name 

RA 

(J2000) 

DEC 

(12000) 

ObsIDs 

Fv 70yum 

[mly] 

Fv 160yum 

[mly] 

SpT 

Cloud 

Ref 

SD 

Multiplicity 

41 

CFHTTAU9 

04 24 26.46 

-h26 49 50.4 

1342227059/60 

10.64 ± 1.09 

< 26 

M6.25 

Tau 

17, 2 



42 

SORI053825 

05 38 25.44 

-02 42 41.3 

1342226721/22 

< 12 

< 22 

M7.0 

o-Ori 

26, 2 



43 

OPHALL1622 

16 22 44.94 

-23 17 13.4 

1342227839/40 

5.75 ± 1.02 

< 160 

M6.0 

Oph 

25,2 

Yl 


44 

2M1160737 

16 07 37.73 

-39 21 38.8 

1342227813/14 

16.15 ± 1.02 

8.54 ± 1.61 

M5.75 

Lup III 

4,2 



45 

DENIS 1603 

16 03 34.71 

-18 29 30.4 

1342227831/32 

3.60 ± 1.01 

< 22 

M5.5 

Upper Sco 

27, 2 



46 

USC0112 

16 00 26.70 

-20 56 31.6 

1342227140/41 

<5 

< 12 

M5.5 

Upper Sco 

28,2 



47 

USC0128 

15 59 11.36 

-23 38 00.2 

1342227134/35 

5.69 ± 1.02 

<25 

M7.0 

Upper Sco 

28,2 

Y2 


48 

USC055 

16 02 45.60 

-23 04 50.0 

1342227136/37 

13.26 ± 1.02 

< 11 

M5.5 

Upper Sco 

28, 2 

Y2 


49 

USD155556 

15 55 56.01 

-20 45 18.7 

1342227138/39 

15.92 ± 1.01 

< 11 

M6.5 

Upper Sco 

27, 2 

Yl, Y2 


50 

USD155601 

15 56 01.04 

-23 38 08.1 

1342227132/33 

7.52 ± 1.01 

<37 

M6.5 

Upper Sco 

27,2 

Y2 


51 

USD 160603 

16 06 03.91 

-20 56 44.4 

1342227142/43 

3.49 ± 1.01 

< 24 

M7.5 

Upper Sco 

27,2 



52 

USD160958 

16 09 58.53 

-23 45 18.6 

1342227146/47 

6.07 ± 1.01 

< 18 

M6.5 

Upper Sco 

27,2 



53 

USD 161005 

16 10 05.42 

-19 19 36.3 

1342227829/30 

5.37 ± 1.01 

8.80 ± 1.61 

M7.0 

Upper Sco 

27,2 



54 

USD161833 

16 18 33.18 

-25 17 50.5 

1342227827/28 

10.42 ± 1.02 

<93 

M6.0 

Upper Sco 

27,2 


binary 

55 

USD161939 

16 19 39.76 

-21 45 35.0 

1342227835/36 

15.88 ± 1.01 

<42 

M7.0 

Upper Sco 

27,2 


binary 


Notes. (1) marks sources without clear excess emission at near- and mid-IR wavelengths. (2) The horizontal line divides the sample into a group of late-type BDs (M8 and later, top) and a group 
of early-type BDs (earlier than M8, bottom). (3) The ObsIDs contain both the scan and cross scan ID numbers of the observations. Taking SSSPM1102 as an example, 1342221849 is the scan ID 
and 1342221850 is the cross-scan ID. All ObsIDs follow the same convention. (4) First reference listed is for spectral type. Further information (e.g., near- and mid-IR photometry) can be found in 
the second reference and references therein. (5) The column “SD” marks the objects for which the re-measured PACS photometry are significantly different from the literature values at 70/jm (Yl) 
and 160;um (Y2 ) respectively , see Se ct. lT^ _ 

Re ferences. (1) IScholz et aU (l2005h : (2) iHarvev et all (l2012ah : (3 ) iGizisI (l2002l) : (4) lAllen et aP (l2007l): (5) jWilking et all (l2005h: (6) lAlves de Oliveira et al.l (l2013l) : (7) IWilking et alj (Il999l) : 
(8) I Alves de Oliveira et al.l ll2010h: (9) iLuh man et alJ ^ 200^: (10 ) iJoergens et al.l fey3h: (1 1) iLi hnianl (120071); (12) ISnezzie t all (l20 03i: (13) ISnezzi et al.l (l201 3h: (14)[Guieu et a l.l (1200' ^; (15 ) 
ICantv et alji2013h: (16) ^Iger et al ](|2014|); (17)lLuhmanl (l20 0^: (18)lL uhman et al.l(l200a): (19) Slesnickrt alj j^Ohh : (20) Briceno et al.lii2002h : (21)lL uhman et al.i( l2009l) : (22) lLuhmanl ( ^2004^ : 
(23) lDamianov et al.l(l2007h : (24) lLuhman et alJ(l2008h : (25) lGullv-Santiago et al.lll201lh : (26) lRigliaco et al.lll201lh : (27) lMarttn et al.l(l2004h : (28) lArdila et all(l200(ih . 
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Table 2. Cloud age and distance assumed in this work. 


Cloud 

Distance 

(pc) 

Age 

(Myr) 

Reference 

Ophiuchus (Oph) 

130 

1-3 

L 2 

Taurus (Tau) 

142.5 

1-3 

3,4 

Chamaeleon I (Cha I) 

162.5 

1-3 

L5 

Chamaeleon II (Cha II) 

178 

1-3 

L6 

Lupus III (Lup III) 

200 

~5 

L 7 

tr Orionis (cr Ori) 

360 

~3 

L 8 

Upper Scorpius (Upper Sco) 

145 

~11 

L9 

TW Hydrae association (TWA) 

55 

8-10 

10, 11 


Referenc e s. (1) iReinurthl dlOOSHl: ( 2) IWilking et all (l2^: (3 ) 
iReinui^ (|2008 ^; (4) iBriceno et alj 1200^: 151 lLuhrMn[ _ ' 
(6) S rezzi et all d200^: ( 7) IComeron et al. | (|2009|); (8) 
yOM ); (9) PecautetalJ (| 2012h : flOl iDucourant et all 
iBarrado Y NavascuesI d2006ll . 


Sherry et al 

20l4l : (111 


3. Observations and data reduction 


All the targets have been observed by the Photodetector 


Array Camera and Spectrograph (PACS; iPoglitsch et al 


on board Herschel with various programs ([Andr e et al 


Harvey et al.ll2012at I Alves de Oliveira et al.ll20l'^ Bulger et alJ 


2010 


201C; 


20141) . The ObsIDs for all the observations are listed in Table [T] 
The majority of the observations were performed in the mini¬ 
scan map mode. The only exceptions are those for GY92 3, 
GY92264, GY92 310 and J125758.7-770120. They were ob¬ 
served in the standard PACS scan map mode. Parallel mode 
observations were not considered. We took the data from the 
Herschel Science Archive at the processing stage of Level 1, 
Version 12.1 of the standard pipeline (SPG = Standard Product 
Generation) and further processed the data by producing maps 
and by conducting photometry. 


3.1. Mapping 

The Level 1 da ta from the archive are fully calibrated 
(iBalog et al.ll2014h . Only the last steps to Level 2 and 2.5 are 
missing, which consist of deglitching, signal drift suppression 
and mapping to convert the detector data into a map that super¬ 
imposes the detector signals for any given position on the sky 
covered during the observation. 

Since we are only interested in the point sources, the high- 
pass hltering algorithm was applied to get rid of the background 
and detector signal drifts. This approach only conserves the 
source flux, if an object mask indicates those parts of the data 
that are not considered as background. For very faint objects like 
BDs, the mask produced by the SPG is insufficient. Therefore, 
we produced new and adequate source masks from Level 2.5 
SPG maps (scan and cross-scan pairs are coadded) based on a 
signal-to-noise criterion. 

After deglitching the data as done by the SPG, we produced 
maps using the photProject task. It converts the flux distribu¬ 
tion as measured by the fixed detector pixel grid into a rebinned 
map, where the flux is redistributed on a freely selectable and 
virtual pixel grid. Since the map pixels are usually smaller, dis¬ 
torted and rotated relative to the detector pixels, the flux distribu¬ 
tion must and can be done with great care. There is a “drizzling” 
algorithm, originally invented for th e HST that allows modif ying 
crucial parameters when doing this (iFruchter & Hookl|2002h . 

The SPG uses a so-called drop size that is relatively large. 
This is a scaling factor of the detector pixel sizes that is ap¬ 


plied before projecting them on the final map grid. This leads 
to a redistribution of the source flux on larger areas of the fi¬ 
nal map. This produces prettier, i.e. smoother images due to the 
convolution of neighbouring detector pixel signals. But this also 
averages out noise, which leads to an underestimate of the real 
detector noise. Reducing the drop size can minimize this effect. 
The images look sharper, but also noisier. However, the result¬ 
ing noise properties are much closer to the reality than what one 
gets with large drop sizes. Therefore, the subsequent photome¬ 
try is also more accurate. We adopted 0.1 as the drop size for all 
maps. The pixel sizes were chosen to be 1" at IQ pm and 2" at 
160jum. 

Finally, we combine the map pairs of a so-called scan and a 
cross scan, corresponding to the Level 2.5 stage of the SPG. 


3.2. Photometry 

Since the sources are very faint and the vicinities around the tar¬ 
gets of interest are full of other and brighter sources, a simple 
aperture photometry procedure is not sufficient. Therefore, we 
applied a HIPE built-in source extraction and photometry tool 
called Sussextractor. It is based on Bayesian statistics and 
delivers the most likely result for the target position and the inte- 
grated source flux by fitt ing 2D Gaussians to the detected sources 
dSavage & 01ivedl2007t) . 

A Gaussian source model is not adequate for the more com¬ 
plicated PSF shape of HerschelfPACS point sources. Only the 
core of the PACS PSF is close to Gaussian, but it also contains 
a trilobal shape at a level of a few percent of the peak inten¬ 
sity superimposed by an extended and slowly declining halo. 
Therefore, a correction factor has to be applied to the photometry 
derived by Sussextractor. 

For this purpose, we reprocessed standard star data of alpha 
Boo (a PACS prime standard) in the same way as the BDs and 
then applied the Sussextractor algorithm and compared the 
results with the ones derived from a standardized aperture pho¬ 
tometry procedure including an aperture correction that was de¬ 
rived during the PACS flux calibration program. The multiplica¬ 
tive correction factors turned out to be 1.513 at IQ pm and 1.427 
at 160 pm. The photometry errors given by the Sussextractor tool 
were corrected in the same way and regarded as the photometric 
errors. Additional systematic uncertainties were not considered. 
For non-detections, upper limits have been obtained by a statis¬ 
tical analysis of a series of apertu re photometry mea surements 
on the background as explained in iBalog et al.l (120141) . This ap¬ 
proach avoids underestimating the noise obtained from the map 
pixel statistics that is affected by correlation artefacts introduced 
by the drizzling algorithm. 

3.3. General trends of the observational data 

The measured flux densities at PACS bands are summarized in 
Table [1] whereas Figure [T] shows a comparison between our re¬ 
sults and previously published photometry. It is clear that most 
detections at 70 pm by other studies are confirmed with our re¬ 
analysis of the data. In particular, our flux measurements of the 
detected sources are in agreement with previous results (see the 
red dots in the upper panel of Figure [T]i. However, we obtained 
less than one half of 160jum detections that are reported be¬ 
fore. This is probably because Herschel 160/tm images of faint 
disks are easily contaminated with the extended emission from 
the background, and the effects of using different data reduction 
algorithms on the results are significant in this case. In Table [1] 
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10 100 1000 
Flux density in this study [mJy] 


Fig. 1. Our HerschelfPA.es flux densities at 70 /mi (upper panel) 
and 160/im (lower panel) as compared to previous results. The 
red dots depict real detections identified by both our and pre¬ 
vious works, whereas green marks consisting of leftward and 
downward arrows symbolize the upper limits. The black squares 
with downward arrows and blue squares with leftward arrows 
show detections by our study and previous works, respectively. 

we mark two kinds of objects for which the re-analyzed PACS 
photometry are significantly different from the literature values: 
(1) the difference in the flux density between our and previous 
results is larger than 3 cr, where cr refers to the RMS level for 
each object given here; (2) firm detections reported before, but 
upper limits deduced by us, and vice versa. Since the constraints 
on the disk structure, especially the structure of the outer and 
cool regions of the d i sk, are sensitive to far -IR data points (e.g., 
iHarvev et al.l2012^llbLlOlofsson et al.l2013h . the discrepancies in 
the PACS photometry between our and previous results will re¬ 
sult in different model explanations of the data. Therefore, our 
homogeneous analysis of the Herschel data and subsequent self- 



70 micron flux density at lOOpc (mJy) 

Fig. 2. The cumulative distribution of the 70/im flux density 
for early-type brown dwarfs (ETBDs, blue) and late-type brown 
dwarfs (LTBDs, red). All flux densities are scaled assuming a 
distance of 100 pc. The median values are indicated as vertical 
dashed lines. 


consistent modeling have the potential to retrieve the hidden sta¬ 
tistical trends of disk properties in the substellar mass regime. 

We generated the cumulative distributions /(> Fy) of PACS 
photometry, using the Kaplan-Meier product-limit estimator to 
properly account for censo red data sets, i.e., upper limits on Fy 
(iFeigelson & Nelsonlll985h . This was done to examine any sta¬ 
tistical difference in the far-IR emission between the two groups 
differing in terms of SpT. All flux densities are scaled by assum¬ 
ing a distance of 100 pc. The results are shown in Figure |2] The 
LTBD sample (red line) has systematically lower flux densities 
than those of the FTBD group (bl ue line), which is consistent 
with the finding bv iLiu et^ (1201 5h that the amount of Herschel 
far-IR emission correlates well with the SpTs of the host BDs. 
From the distribution, the median flux densities are estimated, 
which we show as the vertical dashed lines. The median 70/rm 
flux density of the FTBD group is Fjo - 22 mJy, which is at least 
three times larger than that of the LTBDs, i.e., Fyo < 7 mJy. We 
consider the latter case as an upper limit because more than half 
of the sources are undetected at this wavelength (iMathews et al.l 
l2013h . Due to the same reason, we do not perform the Kaplan- 
Meier product-limit estimation for the 160 /rm data. 

From a theoretical point of view, disks with a larger flar¬ 
ing index intercept a larger portion of the central star’s radi¬ 
ation due to a lar ger flaring angle and cons equently re-emit 
more IR flux (e.g.. Ichiang & Goldreichiri997h . The systemati¬ 
cally lower flux density at far-IR wavelengths of the LTBDs im- 
plies that disks arou nd lower mass stars are generally less flared. 
iBulger et al.l (1201 4ll investigated the effect of individual disk pa¬ 
rameter on the appearance of the SED, and their results indi¬ 
cate that the flaring index and scale height work together to de¬ 
termine the emission level at HerschelfPACS bands. Moreover, 
the observational trends of far-IR flux densities may also help 
us to consider whether there is any evidence for differences in 
disk mass with different stellar masses, since the HerschelfPACS 
photo metry can provide a rough mass estimation of BD disk s 
(e.g., iHarvev et al.l 12012^ ISpezzi et'^l2013L iLiu et al.l 1201 5h . 
Nevertheless, interpreting the observational trends presented 
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here needs detailed SED analysis, in which efforts of reducing 
model degeneracy should be made as far as possible given that 
the data available are limited. 


4. SED Modeling 

All the targets have been observed over a broad range of wave¬ 
lengths. We constructed the broadband SEDs of our sources 
by using our HerschelfPACS flux measurements and adding 
ancillary data at optical and mid-IR wave lengths, from SPSS 
catalo g dAhn et alJ 20121). DENIS survey dPENIS Consortium! 


4.2. Disk model 

Dust distribution: We employed the standard flared disk model 
with well-mixed gas and dust, which has been successfully used 
to explain the observed S EDs of a laree sample of young stel¬ 
lar objects an d BDs (e.g.. Wolf et al.| 2003 t ISauter et al.lj2009l : 
iHarvev et al.ll201^ Lloergens et al.l 2()13t On et alJ 20151) . The 


structure of the dust density is assumed with a Gaussian vertical 
profile 


Cutri et al. 


Pdust “ PO 1 


120051). 2MASS (Cutri et ani2003h . Spitzer and WISE catalogs 


/Ra“ 

[ 1 / z 


2 \h(m)l 


( 1 ) 


120121) . We also took into account the up- and the surface density is described as a power-law function 


(e.g., 

per limits of flux measurements in the (s ub-)millimeter regime 
available for a few of our so urces fe.g.. iMohantvet ani2013t 
iBroekhoven-Eiene et al.l2014l) . We derived the (sub)stellar prop¬ 
erties of each target from thoroughly modeling the photosphere. 
Eor the 46 targets that show IR excess produced by circumstel- 
lar dust, we conducted detailed S ED an a lysis u sing the radiative 
transfer code MC3D developed bv IWolfI (l2003h in order to char¬ 
acterize the structure of their surrounding disks. 


4.1. Determination of the (sub)stellar properties 
We performed a thorough modeling of the photosphere of each 




( 2 ) 


where m is the radial distance from the central star measured in 
the disk midplane, and !i(m) is the scale height of the disk. The 
disk extends from an inner radius Rin to an outer radius Rout- To 
the best of our knowledge, among our sample, there are 5 ob¬ 
jects that have been identified as binary s ystems so far. They are 
2M1207 fa--55 AU.IChpvin et al.ll2()03) . J04221332-H1934392 
(a~l AU.lTodorov et al.ll20l4 . J04414489-H2 301513 (a~15 AU , 


We perrormea a tnorougn moaelmg or me pnotospnere o 
object by fitting the BT-Settl models d Allard et alil2012l) with a 

mi nimization as well as Bayesian statistics (iBavo et ^l2008l 
l2014l) . The BT-Settl models incorporate a sophisticated treat¬ 
ment of photospheric dust, which is likely to affect the atmo¬ 
spheres and correspondingly the radiation field of cool objects 
like BDs in our sample. The broad-band photometry for the pho¬ 
tospheric part of the SED are taken into account in the fitting 
procedure. A similar analysis was done bv iJoergens et alJ (l2013h 
for their case study of OTS44, a BD at the planetary mass bor¬ 
der. 

Eor objects in star-forming regions, we use the standard dis¬ 
tances wit h errors to the clo ud as given by the Handbook of Star 
formation (lReit)urthll2008bll3) . while the kinem atic distances for 
objec ts in loose associations are adopted (e.g.. iDucourant et al.l 
l2014l) . The extinction is treated as a free parameter. Eor objects 
in star-forming regions, we use as starting poi nt upper limits for 
the ex tincti on based on extinction m aps from iKainulainen et alJ 
(l2009l) and ICambresv et al.l (Il997l) . For objects in loose asso¬ 
ciations, we consider an upper limit of 0.1 mag, and for cr Ori, 
we convert the standard E(B-V) value to extinction Ay assuming 
Rv = 3.1 since this is the basic as sumption for the extinction law 
used in VOSA dBavo et al.ll2008l) . 

The derived best-fit (y^) and most probable values for the 
effective temperatures, luminosities, and extinctions are in good 
agreement with each other for ~90% of the objects. For a few 
cases, we only get moderate fitting results: (i) for object 21 
(J04274538-i-2357243),the extinction Ay is loosely constrained, 
we therefore took the value that minimizes the y^. (ii) for ob¬ 
ject 34 (2MJ110703), there is no data point bluer than 2MASS/J 
band. In this case, we fixed Ay to 15 mag which is l ower than 
the value of 16.5 mag estimated in iFurlan et al. ( 20091). bu t is in 
agreement with the one derived bv iLuhman et al. ( 2^008h . (iii) 
object 36 (ISO 252) for which the best and second best fits are 
indistinguishable in terms of;^f^ and the Bayesian analysis shows 
that it is hard to break the degeneracy between Tgff and Ay. We 
took the best-fit parameter set for this object. The derived effec¬ 
tive temperatures, luminosities and extinctions are summarized 
in Table [3] 


iTodorov et al.l20l4 . USD 16 1833 (a~ 134 AU, 


Bouv et al.l2006h 


and USD161939 (a~26 AU. iBouv et al.l 12006 ). where a refers 


to the separation within the system. The disks around individ¬ 
ual components in binary systems are e xpected to have trun¬ 
cation radii of the order of (0.3 - 0.5)a (iPanaloizou & Pringlel 
Il977l) . We adopted 0.5 a as the disk outer radii for 2M1207, 
USD161833 and USD161939. For the close pairs (a< 15 AU, 
J04221332-H1934392 and J04414489-k2301513), dynamical 
simulations of star-disk inter actions suggest that individua l disks 
are unlikely to survive (e.g.. lArtvmowicz & Lubov^l 19941) . Disk 
modeling is complicated in those close multiple systems. For 
simplicity, we assume that the emission is associated with cir- 
cumbinary disks of 100 AU in size. For other objects, we fix 
f^out = lOOAU in the modeling, because the choice of this pa¬ 
rameter value makes essentially no diffe rence to the syntheti c 
SEDs in the simulated wavelength range (iHarvev et al.ll201^ . 
The scale height follows the power law distribution 


! VJ \P 

h(m) = Hwo - , 

^ ^ UOOAU; 


(3) 


with the exponent p characterizing the degree of flaring and 
//lOo representing the scale height at a distance of lOOAU 
from the central star. The indices a, p, and p are codependent 
through p - a - p. We fix p = 1 that is the ty pical value 
found for T Tauri disks in t he sub-millimeter (e.g.. llsella et al.l 
1200^ iGuilioteau et akllMlh . since only spati ally resolved data 
can p lace constraints on this parameter (e.g., iRicci et aH 1201 31 
120141) . 

Dust properties: We assume the dust grains to be a ho¬ 
mogeneous mixture of 75% amorphous silicate and 25% carbon 
with a mean density of p^ rum - 2.5 g cm~^ a nd the complex 


Lgc _ _ 

refract ive i ndices given bv |.Tager et alT d 19941) . iDorschner et al.l 
(1 19951) and Llager et al.l (Il998l) . Porous grains are not considered 
because the fluxes at wavelengths beyond ~2pm are almost 
independent o f the degree of grain po r osity i n low mass disks, 
as shown by iKirchschlager & WolfI (l2014l) . The grain size 
distribution is given by the standard power law d«(a) oc a“^ ^da 
with minimum and maximum grain sizes Omin = 0.1 pm and 


7 
















































































































Liu et al.: Homogeneous analysis of disks around brown dwarfs 


Table 4. Parameter range of the model grids computed in this 
work. 


Parameter 

Min 

Max 

Number of values 

Sampling 

^in sub] 

1 

100 

6 

logarithmic 

P 

1.0 

1.3 

13 

linear 

Eioo [AU] 

2 

20 

10 

linear 

log Mdisk [Mq] 

-6 

-2.5 

8 

logarithmic 

in 

15 

90 

6 

linear 


Qmax = lOOyum, respectively. The choice of the minimum 
value for the grain size, amin, ensures that its exact value has 
a negligible impact on the synthetic SEDs. Since there is no 
information about the maximum grain sizes of our target disks, 
as provided e.g., by the (sub-)millimeter spectral index, we 
adopt amax = lOOjum. The HerschelfPACS far-IR observations 
are sensitive to dust grains with this assumed sizes. Strong 
grain growt h up to mil l imeter sizes a s detected in some BD 
disks (e.g., RiccigLalJ 1201 2L 1201 3t iBroekhoven-Fieneet alJ 
l2014t iRicci et ^ 2014 ) would remain undetected in our data 
and could affect the disk mass. Our prescript i on for the dust 
properties is identical to those used in lLiu et SI (l2015h . 


Heating sources: The disk is assumed to be passively 
heated by stellar irrad iation (e.g., IChiang & Goldrei^ Il997t 
iDullemond et alJ 1200ih . Other heating sources, such as the 
viscous accretion, are not taken into account as this would 
only introduce further free parameters without any qualitative 
constrai nts from the modeli ng. We use the BT-Settl atmosphere 
models (lAllard et al.l 12012h as incident flux with parameters 
as constrained in Sect. 14.11 The radiative transfer problem is 
solved self-consistently considering 100 wavelengths, which are 
logarithmically distributed in the range of [0.05 /tm, 2000 gLm\. 

4.3. Fitting method 

The task of SED fitting was conducted with a grid of pre¬ 
calculated models for each target. The model grids can give 
us an overview of the fitting quality in different parameter do¬ 
mains and clues to improve the quality of the fit with further 
attempts if necessary. Table |4] summarizes the explored ranges 
of each parameter considered here. The grid adds up to a total 
number of 37 440 models for each source. We tried to achieve 
a better explanation on the observations on one hand by reduc¬ 
ing the dimensionality of the parameter space as far as possi¬ 
ble and on the other hand by using a denser grid of models 
as compar ed to previous studie s on HerschelfPACS data of BD 


disks (e. ff.. IHarvev et al. 2012allbl: 

Alves de Oliveira et al.ll2013l 

lOlofsson et al.ll2013l:ISDezzi et al.l 

2013). For obiects. for which 


we could not obtain good results from the grid, simulated an¬ 
nealing approach was additionally invoked to improve the solu¬ 
tion by taking the best model in the grid as the starting point 
of the Markov chain. Generally, ~ 1 000 steps for t he local fit 
with simulated annealing are enough (iLiu et al.l2013h . This kind 
of methodology makes use of the advantages of both database 
method and simulated annealing algorithm and has already been 
demonstrated to be su ccessful for SED analysis (iLiu et al.l2012t 
iMadlener et al.l[^12h . 

The fitting results are displayed in Figure [3 The best-fit 
models are indicated as black solid lines, whereas the dashed 
lines represent the photospheric emission level. The correspond¬ 


ing disk parameter sets are listed in Table [3 The SpitzerflRS 
spectra, when available, are taken into account by fitting the 
subjacent continuum. Since our goal is to characterize the over¬ 
all structure of the disks, we therefore did not attempt to repro¬ 
duce the exact shape of the silicate feature that is mainly related 
to silicate mineralogy of the dust in the di sk atmosphere (e.g., 
iBouwman et al.ll3i)08t lOlofsson et al.ll201()h . We emphasize that 
the best fits presented here cannot be considered as a unique so¬ 
lution due to the degeneracy of SED models. Despite this, pre¬ 
vious studies have shown that modeling SEDs with broad wave¬ 
length coverage c an constrain the mass and geometry of disks 
around BDs (e.g. , iHarvev et al.l l2012albt lOlofsson et al.l 120131 : 
ISnezzi et al.ll2013h . In particular, the Bayesian inference can an¬ 
alyze the potential correlations and i nterplay between different 
parameters in a statistical way (e.g.. iPinte et alj|2008l) . We de¬ 
rived the Bayesian probability distributions of the disk param¬ 
eters, from which we calculated the most probable values and 
estimated the validity ranges for each parameter, corr esponding 
to regions where P > 0.5 x EMax (e.g.. lLiu et al.ll20f3i . The gray 
lines in Figure[3denote all the models that are within the validity 
ranges. 

5. Discussion 

We re-analyzed the HerschelfPACS data of a sample of 55 BDs 
and very low mass stars in a homogeneous way and constructed 
their observed SEDs by complementing our Herschel photom¬ 
etry with previous flux measurements available at wavelengths 
ranging from optical to (sub-)millimeter. For the 46 targets that 
show IR excess emission, we characterized the properties of 
their surrounding disks using radiative transfer technique un¬ 
der identical assumptions for the modeling setup. Moreover, we 
evaluated the constraints on different disk parameters through 
Bayesian analysis. 

Despite the good wavelength coverage of observations and 
efforts to reduce the degree of freedom, the degeneracy of SED 
models still exists. Therefore, we statistically analyzed the most 
probable results from Bayesian statistics, other than the best-fit 
models, because the former case is likely more representative 
of the overall fitting quality. The most probable parameter sets 
for each target are identified as the peak of Bayesian probability 
distribution. The derived values are summarized in Table[3 

5.1. Overview of modeiing resuits 

We characterized the studied disks by determining their inner 
radii Ein, flaring indices /?, scale heights //loo, disk masses Mdisk 
and inclinations i. The results are listed in Table [3 We give an 
overview of the statistical behaviour of these properties in the 
following. 

The disk inner radius can be well constrained for all objects 
because the observed SEDs are well sampled in the near- and 
mid-lR domains. The inner radii of the disks are likely close 
to the dust sublimation radii rsub, which are of the order of 
[0.005-0.05 AU] from the faintest to the brightest targets in our 
sample. We found four candidates of BD transition disks with 
Ein > 20rsiib, i.e., source 22 (J04302365-b2359129), source 24 
(J04355143-b2249119), source 25 (J04361030H-2159364) and 
source 29 (J04221332-bl934392). Current data indicate no de¬ 
tectable excess shortward of 12/im accompanying with a steep 
slope in the mid-lR of these four objects. However, none of them 
is detected at HerschelfPACS bands. Follow-up observations, in 
particular spatially resolved millimete r images, are required to 
clarify their evolutionary stages (e.g.. IWilliams & Ciezal 120111 : 
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Table 3. Photospheric and best-fit/most-probable disk parameters from this work. 


ID 

Tcff 

L. 

Av 

Rm 

Afdisk (B) 

Afdisk (P) 

P{B) 

PiP) 

o 

o 

Hioo(P) 

i 


[K] 

[10-^ Lo] 

[mag] 

[AU] 

1 

o 

1 

o 



[AU] 

[AU] 

r] 

1 

2400 

0.355 ±0.106 

0.1 ±0.01 

0.004 

1.0 

1.0 

1.10 

1.10 

20 

20 

60 

2 

2500 

2.367 ± 0.336 

0.1 ±0.01 

0.025 

0.1 

0.1 

1.175 

1.175 

14 

14 

45 

3 

2100 

6.139 ± 1.852 

0.26 ± 0.26 

0.038 

3.16 

1.0 

1.10 

1.10 

16 

16 

15 

4 

2800 

64.2 ± 10.49 

3.6 ±0.15 

0.041 

3.16 

3.16 

1.10 

1.10 

16 

16 

30 

5 

2800 

23.18 ±4.020 

2.28 ± 0.18 

0.025 

10.0 

10.0 

1.10 

1.125 

20 

18 

30 

6 

2800 

81.15 ± 12.94 

7.2 ± 0.25 

0.317 

10.0 

10.0 

1.125 

1.125 

14 

14 

60 

7 

2600 

379.1 ±58.95 

3.94 ± 0.98 

0.015 

0.1 

0.32 

1.15 

1.15 

20 

20 

15 

8 

1700 

2.4 ± 0.54 

2.6 ± 0.20 

0.023 

3.25 

3.16 

1.31 

1.30 

15.6 

16 

60 

9 

2800 

5.459 ± 0.352 

3.0 ±0.19 

0.038 

3.16 

3.16 

1.20 

1.175 

8 

10 

15 

10 

2500 

6.242 ± 0.325 

4.8 ± 0.22 

0.005 

0.1 

0.1 

1.025 

1.05 

18 

18 

60 

11 

2600 

3.54 ±0.192 

1.35 ±0.1 

0.005 

1.0 

0.32 

1.10 

1.10 

8 

10 

30 

12 

2600 

5.603 ± 0.376 

1.41 ± 0.07 

0.038 

1.0 

1.0 

1.10 

1.10 

8 

8 

60 

13 

2200 

1.221 ± 0.080 

0.2 ± 0.03 

- 

- 

- 

- 

- 

- 

- 

- 

14 

2400 

3.22 ±0.152 

1.27 ± 0.07 

- 

- 

- 

- 

- 

- 

- 

- 

15 

2500 

5.591 ± 0.237 

1.89 ±0.1 

- 

- 

- 

- 

- 

- 

- 

- 

16 

2600 

4.078 ± 0.177 

3.91 ± 0.23 

0.010 

100 

100 

1.10 

1.075 

20 

20 

45 

17 

2700 

10.43 ± 0.429 

1.71 ± 0.1 

0.051 

3.16 

1.0 

1.10 

1.125 

8 

10 

45 

18 

1800 

3.045 ±0.185 

1.04 ±0.1 

- 

- 

- 

- 

- 

- 

- 

- 

19 

2500 

3.749 ±0.159 

3.45 ± 0.18 

- 

- 

- 

- 

- 

- 

- 

- 

20 

2400 

4.044 ±0.181 

0.42 ± 0.42 

0.008 

100 

100 

1.137 

1.125 

21 

20 

45 

21 

2200 

2.814 ± 0.165 

0.57 ± 0.08 

- 

- 

- 

- 

- 

- 

- 

- 

22 

2200 

2.66 ±0.165 

0.39 ±0.13 

0.4 

3.16 

1.0 

1.0 

1.0 

11 

12 

30 

23 

2700 

18.74 ± 0.799 

2.42 ±0.17 

- 

- 

- 

- 

- 

- 

- 

- 

24 

2400 

2.517 ± 0.139 

2.14 ±0.16 

0.4 

31.6 

31.6 

1.05 

1.05 

20 

18 

60 

25 

2000 

2.8 ± 0.140 

0.0 ± 0.04 

0.159 

3.16 

1.0 

1.05 

1.075 

8 

6 

45 

26 

2400 

3.564 ±0.171 

0.4 ±0.1 

- 

- 

- 

- 

- 

- 

- 

- 

27 

2300 

2.862 ± 0.137 

0.9 ±0.1 

- 

- 

- 

- 

- 

- 

- 

- 

28 

3000 

1.459 ± 0.076 

2.89 ± 0.16 

0.008 

3.16 

0.32 

1.0 

1.025 

11 

12 

30 

29 

2200 

20.46 ± 0.884 

0.56 ± 0.04 

0.9 

1.0 

0.1 

1.0 

1.0 

2 

2 

60 

30 

2300 

5.669 ± 0.237 

1.43 ± 0.07 

0.013 

3.16 

3.162 

1.05 

1.025 

11 

8 

45 

31 

2800 

10.15 ± 0.612 

0.95 ± 0.32 

0.044 

0.32 

0.32 

1.20 

1.20 

18 

20 

45 

32 

2800 

47.78 ± 1.954 

4.7 ± 0.36 

0.0377 

10.0 

10.0 

1.05 

1.075 

8 

10 

45 

33 

2800 

82.84 ± 3.77 

2.8 ± 0.2 

0.126 

10.0 

10.0 

1.125 

1.125 

4 

4 

45 

34 

2600 

27.39 ± 0.995 

15.0 

0.051 

316 

316 

1.175 

1.175 

20 

20 

45 

35 

2800 

51.18 ± 2.105 

5.42 ± 0.42 

0.041 

3.16 

3.16 

1.20 

1.20 

10 

12 

45 

36 

1900 

17.52 ± 0.933 

2.31 ±0.33 

0.051 

10.0 

3.16 

1.025 

1.05 

12 

12 

45 

37 

3000 

64.0 ± 2.775 

3.58 ± 0.36 

0.017 

3.16 

3.16 

1.25 

1.25 

18 

18 

45 

38 

2200 

30.73 ± 6.521 

7.74 ± 0.48 

0.025 

0.32 

1.0 

1.10 

1.125 

18 

18 

60 

39 

2400 

26.27 ± 5.909 

2.3 ± 0.21 

0.063 

1.0 

1.0 

1.15 

1.15 

10 

10 

45 

40 

3000 

43.82 ± 1.729 

4.3 ± 0.29 

0.015 

0.1 

0.1 

1.05 

1.025 

4 

4 

30 

41 

2800 

23.81 ± 1.39 

1.16 ±0.1 

0.063 

0.32 

0.32 

1.175 

1.175 

18 

16 

60 

42 

2800 

28.93 ± 5.848 

0.04 ± 0.01 

0.025 

10.0 

10.0 

1.075 

1.075 

18 

14 

60 

43 

2100 

12.48 ± 2.211 

1.93 ±0.19 

0.044 

0.32 

0.32 

1.225 

1.225 

18 

18 

30 

44 

2400 

16.49 ± 5.241 

1.02 ±0.26 

0.018 

1.0 

1.0 

1.225 

1.225 

20 

18 

30 

45 

2800 

34.59 ± 7.973 

1.0 ±0.1 

0.013 

0.1 

0.1 

1.275 

1.275 

10 

10 

30 

46 

3000 

13.02 ± 3.103 

0.6 ±0.1 

0.021 

0.32 

0.32 

1.175 

1.175 

12 

12 

60 

47 

2500 

6.242 ± 1.354 

1.3 ±0.1 

0.032 

3.16 

10.0 

1.175 

1.175 

16 

16 

30 

48 

2900 

30.6 ±7.081 

0.3 ±0.1 

0.012 

1.0 

1.0 

1.15 

1.15 

14 

14 

45 

49 

3400 

27.66 ± 6.102 

3.38 ± 0.23 

0.010 

1.0 

1.0 

1.20 

1.20 

16 

16 

45 

50 

2800 

10.34 ± 2.422 

1.3 ±0.1 

0.007 

31.6 

31.6 

1.125 

1.15 

18 

18 

30 

51 

2800 

14.88 ± 3.411 

1.5 ±0.1 

0.051 

1.0 

1.0 

1.10 

1.10 

6 

6 

60 

52 

2600 

26.31 ± 6.167 

0.6 ±0.1 

0.01 

0.32 

0.32 

1.20 

1.225 

10 

10 

45 

53 

2600 

6.486 ± 1.574 

0.7 ±0.1 

0.006 

10.0 

10.0 

1.10 

1.10 

20 

20 

60 

54 

3000 

42.35 ± 9.502 

1.7 ± 0.1 

0.014 

0.3 

0.32 

1.20 

1.20 

8 

8 

30 

55 

3000 

21.45 ±4.865 

1.5 ±0.1 

0.025 

0.15 

0.32 

1.175 

1.175 

20 

18 

45 


Notes. (1) We list the disk parameters for all objects in our sample with detected excess emission at IR wavelengths. (2) The error in Tcs is 100 K 
for all cases except for object 8 (i.e., OTS 44) with a larger value of 140 K. (3) We fixed Av/ = 15 mag for object 34, see an explanation in Sect. l4.ll 
(4) The total disk mass is calculated from the dust mass assuming a gas-to-dust mass ratio of 100. (5) The tags B and P quoted in brackets 
for Mdisk, P and Hioa denote the best-fit and the most probable values, respectively. The values of the disk inner radius R\^ and the inclination i are 
for the best-fit models. 
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Fig. 3. Spectral energy distributions of the target disks. The dots depict the photometry at various wavelengths. Spifzer/IRS spectra 
are indicated as red lines. The upside down triangles show the 3 cr upper limits of the flux density. The best-fit models are indicated 
as black solid lines, whereas the dashed lines represent the photospheric emission levels. The gray lines denote all probable model 
fits that are derived by Bayesian analysis, showing the strength of constraints of the best fit models for each object, see Sect. 14.31 


[Andrews et alJ 1201 ih . The inclination of the disk cannot be 
constrained very well with the available SED data as shown 
by the flat probability distribution of this parameter nearly for 
all targets. Accurate determination of disk inclinations requires 
spatially resolved obs ervations (e.g.. lChiang & Goldreichlfl999t 
iRobitaille et al.ll200^ . 

The geometric structure of the disk is described by the flar¬ 
ing index (3 and scale height //mo according to the model ansatz, 
see Sect. 14.21 Our detailed SED analysis places tight constraints 


on the disk geometry, because the best-fit p and //mo are in good 
agreement with the most probable values (see Table [3]) and the 
probability distribution functions are clearly non-flat for most 
objects (see the Appendix). We observed a preferential values 
of p in the range 1.05 - 1.20, while the majority of the target 
disks are best modeled with // mn o f the order of 5 - 20 AU. As 
shown bv iHarvev et al.l (12012ah and lLiu et al.l (12015h . most parts 
(except for the inner region) of typical BD disks produce op¬ 
tically thin emission at far-IR wavelengths, demonstrating the 
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Fig. 3. Continued. 


applicability of PACS photometry as a diagnosis for the disk 
mass of BDs. Objects with detections at 70 and/or 160m al¬ 
low reasonable determination of their disk masses as demon¬ 
strated by clear peaks in the Bayesian probability distributions 
of Mdisk in these cases. For the remaining targets, the disk mass 
is not well constrained, typically showing two or three indis¬ 
tinguishable peaks in the probability distributions. Overall, the 
likely disk masses feature a wide range of [10“®, 10 “"^ Mq] with 
a median value of the order of 10 “^ Mq, where we assume 
the dust properties as described in Sect. 14.21 and a gas-to-dust 
mass ratio of 100. Our statistical results on the disk param¬ 
eters of substellar hosts are widely in agreement with previ¬ 


ous findings from case studies and surveys, however, for indi¬ 
vidual sources presented here our measured PACS flux densi¬ 
ties and correspondingly t he model solutions d iff er from those 


or previous worKS le.s., j-iarvev ei ai.i zuizayoc iKiaz (X uizisi 
2012; lAlves de Oliveira et al. 2013; Joereens et al.ll2012, 2013t 

Snezzi et al. 2013). A 

'ew BDs detected at (sub)millimeter to 
xplanations also support our results on 
operties especiallv the mass of BD disks 

date and their model e 
the tvpical structural pr 

("e.g.. Klein et al.ll2003l 

Scholz et al.ll200d Mohantv et al.ll2013t 

Ricci et alJl2012L 

2013 

2014l;lBroekhoven-Fiene et alj|2014ll. 
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Fig. 3. Continued. 


5.2. Comparison with disks around young steiiar objects 

Comparing the derived disk properties with those of their 
higher mass counterparts, such as T Tauri disks, can provide 
important insights into understanding the formation of BDs 
and very low mass s tars. Several attemp ts have been mad e 
in this direction (e.g., iPascucci et al.ll200^ ISzucs et ^l2010h . 
Plenty of multi-wavelength data including spatially resolved 
ones exist for young stellar objects, leading to robust con¬ 
straints on the structure and mass of their disks. In partic¬ 
ular, coherent multi-wavelength modeling suggests a typical 
flaring index /3 ~ 1.25 and a favorable range 5 - 20 AU of 
H\CIO for T Tauri disks (e.g., Wolf_eta^ 2003 Walkeret^ 


200llSauter et al.l2009l:lMadlener et alj2012HGrafe et alJ2013t 


the ranse of TIO 10 ^ fe.ff. lAndrews & WilliamsI 12007 

Lee et al.ll2011’ 

Sicilia-Aguilar et al.l 201 ll: lAndrews et al.l2013 

Carpenter et all 

201411. 


BDs and very low mass stars are generally flatter and orders 
of magnitude less massive than T Tauri disks. Although simi¬ 
lar conclusions were drawn by previous studies (IPascucci et al.l 
120091 : ISzucs et ^ 120101 : lOlofsson et al.l 12013h . the evidence 
we found should be more significant because our results are 
based for the first time on a homogenous analysis all the way 
from Herschel data reduction and flux density measurement to 
SED modeling, and to interpretation. Theoretical models pre¬ 
dict that disks around cooler stars should be more extended in 
the vertical direction (Walker et al. 200 4). However, our re¬ 
sults, together with previous studies (e.g., Harvey et ^l2012^ 


Garufl et^l2014l) . Based on (sub-)millimeter measurements of 


a large sample of young stellar objects in various star forma- 


We found significant evidence for different properties of 
disks between sun-like and cool stars, showing that disks around 


lAlves de Oliveira et al.l201^IOlofsson et ai.l2013[) . indicate that 

the disk scale height is independent of the host stellar mass. Dust 
growth and settling are the two pivotal physical processes simul¬ 
taneously shaping the disk structure, for instance reducing the 
disk scale height (ID’ Alessio et al.ll2006l) . Mid-IR observations 
have shown that these p rocesses probably occur in early stages o f 
the disk evolution (e.g. jFurlan et al.l2006l:lMcClure et al.l2010ll . 
Given the broad age spread, our targets may be at different evo¬ 
lutionary stages during which the degrees of dust growth and 
settling are not at the same level. This will probably weaken any 
relation, if present, between the scale height and the (sub-)stellar 
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mass. Nevertheless, analyse s of coeval s amples, for example in 
Chal (lOlofsson et al.ll2()l^ and p Oph (lAlves de Oliveira et al.l 
1201 3h . also indicate that disk scale heights are comparable 
in both sun-like stars and BDs. Future high-resolution (sub- 
)millimeter observations are indispensable to quantify the degree 
of dust settling and better constrain th e disk scale height (e.g., 
ISauter et'ani2009^ iBoehler et al.ll2013ll . 


5.3. Dependence of disk parameters on the mass of 
substellar hosts 


There is growing observational evidence that disk evolution de¬ 
pends on the mass of the central object. For example, the pri¬ 
mordial dust disks are believed to live longer for cool stars, 
because the observed disk frequency of cool s tars and BDs is 
higher in comparison with sun-like stys (e.g.. [Carpen ter et al.l 
2006t LScholz et al.ll2007t fRiaz & GizisIbOOSHBavo et al.ll2012t 


Ribas et al.l 1201 5h ~ The observed trend of small disk accre¬ 

tion rates of BDs indicates that the scaling between the accre¬ 
tion rate and stellar mass (i.e., Macc deduced mainly 

from low- and inter mediate-mass stars extends to the substellar 
mass spectrum (e.g., |Muzerolle e t al.l2005l:lMohantv et al.l2d0^ 


iHerczeg et al.ll20d^ iBavo et ^ 2012h . Analysis of Spitzer/IRS 

spectra demonstrates that the dust processing, such as grain 
growth and crystallization, in the mid-IR emitting regions of 
BDs and very low mass stars appears to be more advanced 
than in disks a r ound sun-like stars with simila r age (e.g., 
Apai et aP 2005t Pascucci et al.l 120091: iRiazI 1200^ . Moreover, 


Pascucci et al.l (1201 3h also found stellar-mass-dependencies in 


the atomic and molecular content of disk atmospheres. 

Searching for evidence of stellar-mass-dependent disk prop¬ 
erties is very common and interesting, be cause it has impor tant 
implications for planet formation theories. fSzucs et di¬ 

vided -200 young stellar objects and BDs in the Chamaeleon I 
star-forming region into two groups with SpTs earlier or later 
than ~M4.5. Through modeling the median SEDs (longward of 
MIPS 24 pm) of these two groups, they found that the disk is on 
average flatter in the lower stellar mass group than in the group 
with higher stellar masses. Based on Spitzer data, s everal other 
works hinted at similar results (e.g., iRiaz et alj2012h . As laid out 
in Sect. 15.21 our results derived by including the Herschel far-lR 
observations evince the same tendency. Note that this is a sta¬ 
tistical behaviour shown from the ensemble of the 46 modeled 
very low-mass stars and BDs. 

As a further step towards examining whether the trends of 
disk properties extend to the low-end of the substellar mass 
regime, we performed a detailed comparison of the fitted disk 
geometry parameters between the ETBD and LTBD groups as 
introduced in Sect. |2] The results are displayed in Eigure ID in 
which the upper panel shows a clear decrease of disk flaring from 
early to late type BDs, with a typical dispersion of the order of 
-0.1. The histogram inside the figure confirms the relationship 
between the flaring index p and SpT. In particular, the median 
value of y6 is 1.10 in the LTBD group, while the ETBD group 
features a higher median value of 1.175. The mean value of (i is 
also found to be smaller in the LTBD group, i.e., jSmean = 1.10, 
as compared to /Jmean = 116 for the ETBD group. The differ¬ 
ence in the median/mean disk flaring between the two groups is 
Ay6 - 0.07, which is about three times of the step size for this di¬ 
mension used to build the model grid. The lower fi of cooler BDs 
suggests that disk flaring is very sensitive to far-IR photometry, 
because the LTBD group is generally fainter at IQ pm, see Sect. 
13.31 OTS44, an M9.5 BD with an estimated mass of -12Mjup 
in the Cha I star-forming region (iBonnefov et al.l2014l) . is an ob- 


1.3 


1.2 

1.1 


1.0 


M5 M6 M7 M8 M9 LO 


25 


20 


- 15 

O 

O 

I 10 


5 


0 

M5 M6 M7 M8 M9 LO 
Spectrol Type 

Fig. 4. Elaring index /? (top) and scale height H\qq (bottom) as a 
function of SpT. We found that (3 decreases clearly from early- 
type BDs to late-type BDs (top), while //loo appears independent 
of the spectral type in the BD regime. The inserted diagrams 
in each panel show the histograms of these two parameters cal¬ 
culated for targets that belong to the ETBD and LTBD groups 
respectively. The vertical dashed lines symbolize the criterion 
(i.e., M8) used to divide the targets, see Sect. |2] Blue dots rep¬ 
resent objects located in clusters with ages <5 Myr, whereas the 
red ones refer to older (>5 Myr) targets. The blue square in the 
upper panel points to the obvious outlier target (OTS 44) of the 
observed relation. Eor better representation, slight changes in ei¬ 
ther the SpT or disk parameters are made for objects with a pair 
of identical values. 


Spectral Type 
Number 

87654321012345 

T-1-1-1-1-1-1-1-1-1-1-1-1-T 


I 

I 

I 



I 

I.I.I.I 


Number 

7654321012345 



vious outlier of the observed relation, see the blue square in the 
upper panel of Eigure|4] The significant IR excess together with 
a rising SED fro m 24 to 70 um require a hig hly flared disk to re¬ 
produce the data (iJoergens et al.l201^l2015i this work). We note 
that our re-measurement of the HerschelfPACS flux densities of 
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Number 

109 8765432101234567 



Spectrol Type 



1 10 100 1000 
\ (jLim) 





Fig. 5. Same as in FigureH) but for the disk mass Mdisk- 


OTS44 l eads to a slightly sm aller disk mass of 3.2 x 10 
(~10Mgi: lJoergens et alj 2 0151 this work) co mpared to a previ¬ 


ous estimate (~3 0Mm; Joergens et al. 2()T^ based on the pho 


tometry given by iHarvev et al.l (l2012a ). VLT/SINFONI spectra 
reve aled active mass acc retion for this extremely low-mass ob¬ 
ject (l.loergens et al.l201^ . suggesting that OTS 44’s disk is prob¬ 
ably at its early stage of evolution, during which dust settling 
(a key role of reducing /3) has not yet proceeded very far. The 
observed SpT - p dependency is unlikely a bias of the mixed 
ages and evolutionary stages of the sample, since both the young 
(< 5 Myr, blue points) and old (> 5 Myr, red points) subgroups 
display the same tendency, see Figure 01 

As shown in the lower panel of Figure 0] there is no clear 
trend visible in the distribution of scale height. The median or 
mean values of //mo are close to each other for both the ETBD 
and LTBD groups, i.e., //ioo~14 AU. This means that the inde¬ 
pendency of //loo on SpT is also found in the low-mass substel- 
lar regime. The disk scale height varies like //oc ((r^j/M*)®-^ if 
the hydrostatic equilibrium between the dust and gas phases is 
assu med, where (Td) is the density weighted mean dust temper¬ 
ature. |Andrew£erai](l20T3l) found that {Tj) scales with the host 
stellar luminosity basically according to {Td) ~ 25(L*/Lo)*^^ K. 
The stellar luminosities L* oc in their sample that covers 

the substellar regime. This suggests that the correlation between 
the scale height and M* may be intrinsically not tight. The scal¬ 
ing between the scale height and M* highly depends on the de¬ 
tailed coupling (i.e., dust settling) between the gas and small 
dust grains, which is quite complicated and has to be investi¬ 
gated withJuturehigl>resolution_(subAmillimeter observations 
(e.g.. [^uter et al.ll2009HBoehler et al]l2013h . 

The median disk masses for both the ETBD and LTBD 
groups are Mdisk = Mp, a value es sentially consistent with 
the one given bv iHarvev et al.l (l2012ah . As shown in Eigure|5] 
we did not find an obvious correlation between the disk mass 
and SpT (as a proxy for the mass of the central object) in the 
BD mass regime, which seems to deviate from the speculation 
from the tight correl ation between th e PACS 70jum flux density 
and SpT shown by iLiu et alJ (l2015h . This is probably due to 
one or more of the following reasons: (1) the constraints on the 


Fig. 6. Left panels: Median SEDs of Class II ETBDs (upper) and 
LTBDs (lower) in the Taurus molecular cloud. The dots depict 
the median flux densities at various bands. The upside down tri¬ 
angles show that the median flux densities at the corresponding 
wavelengths are treated as upper limits. The best fit models are 
indicated as solid lines, whereas the dashed lines mark the pho- 
tospheric emission levels. In the fitting procedure, the disk scale 
height //loo and inclination i are fixed to 14 AU and 45° respec¬ 
tively. The parameter set of the best model in the upper panel: 
Rin - 0.09 AU, Mdisk = 3 X 10“^ Mq, f — 1.16. The parame¬ 
ter values of the best model in the lower panel: - 0.03 AU, 

Mdisk = 5 X 10“® Mq, P - 1.14. Right panels: 2D contour plots of 
the function projected over disk mass Mdisk and flaring index 
p. Contours with different shades of red are drawn at the 68%, 
90%, 95% 2D confidence intervals. The blue stars represent the 
best-fit disk models. The hashed regions correspond to models 
that overpredict the upper limits of the median flux density. 


disk mass are mostly from the IQ pm photometry in our study. 
However, only 52% (i.e., 11/21) of the modeled LTBD objects 
are detected at IQ pm, whereas the detection rate is 88% (i.e., 
22/25) for the ETBD group, see Table [1] The (most probable) 
disk masses of the undetected sources are quite uncertain, po¬ 
tentially preventing the underlying correlation, if present, from 
showing up. (2) The conversion from PACS70/im flux den¬ 
sity to disk mass is not straightforward. On one hand the tar- 
get disks are not com pletely optically thin at PACS wavelengths 
dHarvev et al.ll201^ . and on the other hand IQ pm is not in the 
clear Rayleigh-Jeans tail of the SED. (3) Our correlation anal¬ 
ysis is based on a sample of targets that span a wide range of 
ages. Therefore, the targets are probably at different evolutionary 
stages, which may weaken any relation between the (sub)stellar 
and disk masses. Beside the issue of incompleteness of our sam¬ 
ple selection in each cloud, the detection bias can also affect the 
correlation because more distant clouds have a worse sensitivity 
to disk mass compared to closer siblings. 

In order to further investigate the Mdisk - SpT relation in the 
substellar mass regime, we compile the median SEDs of Class II 
Taurus sources with SpTs lying in the same range defined in 
our sample selection. Our compilati o n is e ntirely based on the 
catalog as reported by iBulger et al.l (l2014h . We only consider 
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sources that have flux density measurements including upper 
limits available at all the 2MASS J/H/K bands, four IRAC and 
WISE bands, MIPS 24/rm, and PACS 70 and 160/tm. With this 
criterion, there are in total 32 targets entering in the analysis. 25 
sources are ETBDs, whereas 7 targets are LTBDs. The Kaplan- 
Meier product-limit estimator is used to properly account for 
upper limits, if encountered, in the estimation of median flux 
density . All the PACS flux densities are taken from iBulger et al.l 
(12014l) for homogeneity, although we re-performed photometry 
for some Taurus targets. The median SEDs of Taurus ETBDs 
and LTBDs are shown as blue dots in Eigure | 6 l An important 
result is that the median E 70 = 36 mJy of Taurus ETBD group 
is at least twelve times larger than that of the Taurus LTBDs, 
i.e., E 70 < 3 mJy. The difference in 70 pm flux density between 
early- and late-type targets becomes significantly larger in co¬ 
eval planet-forming disks as compared to the result from our 
sample of objects with different ages, see Sect. 13.31 We fit the 
median SEDs by using the disk and dust models as described 
in Sect. 14.21 As we discussed above, the disk scale height Hiqq 
is independent on SpT. Hence, we fixed Hiqq to 14 AU, i.e., the 
aforementioned median value of Hiqq of our sample. The disk 
inclination i was set to 45° because SED modeling cannot pro¬ 
vide tight constraints on this parameter. These assumptions were 
made in order to further reduce the model degeneracy, and there¬ 
fore to focus on the parameter study between (3, Mdisk and SpT. 
Eigure 0 shows the fitting results and 2D contour plots of the 
function projected over Mdisk and (3. The flaring index [3 is 
not well constrained in Taurus ETBD group, with probable val¬ 
ues extending to a lower level than the case for Taurus ETBDs. 
This may be an indicator that disks around lower mass objects 
are indeed less flared in general. The best fit value together with 
the confidence intervals clearly indicate that Taurus ETBDs are 
very likely more massive than Taurus LTBDs, which appears 
discrepa nt with the theoretical predicti ons of disk fragmentation 
models (IStamatellos & Herczedl^OL^ . However, further efforts 
are highly desired to clarify this issue given the fact that the disk 
mass presented here is roughly estimated with far-IR data. With 
its unprecedented performance, ALMA can observe very faint 
disks at longer wavelength and provide robust measurement of 
disk mass (e.g.. lRicci et al.l2014li . so as to testify whether the BD 
disks obey the same Mdisk-SpT scaling law to the corresponding 
relations established for higher mass stars. 

6. Summary 

We re-processed the HerschelfPACS data and measured the 70 
and 160/rm flux density for a sample of 55 BDs and very low 
mass stars with SpT ranging from M5.5 to LO. Supplemented 
with previous observations at shorter and (sub-)millimeter wave¬ 
lengths, we constructed their broadband SEDs with extended 
wavelength coverage, providing a valuable opportunity to char¬ 
acterize the disk properties in the substellar mass regime. We di¬ 
vided our sample into two groups, i.e. an ETBD and ETBD group 
that consists of targets with SpT earlier or later than M 8 . We 
found significant difference in the cumulative distributions of the 
flux densities between these two groups. In particular, the ETBD 
group features a median 70/rm flux density at lOOpc that is at 
least three times higher than that of the other group. We studied 
the dependence of disk properties on the SpT (as a proxy for the 
mass) of the central object down to the low-mass BD regime and 
even entering the planetary regime. This knowledge is important 
for us to understand the formation mechanism of BDs and to 
improve the planet formation models. Unlike previous studies, 
our analysis is homogeneous all the way from HerschelfPACS 


data reduction and flux density measurement to SED modeling 
within the substellar regime. 

Eor the 46 objects that show IR excess above the photo- 
spheric emission, we performed detailed SED analysis using 
self-consistent radiative transfer models. We introduced as few 
free parameters as possible in the fitting procedure in order to 
reduce the model degeneracy. The statistics based on the entire 
sample shows that the disk flaring of BDs and very low mass 
stars is indeed smaller than that of their higher mass counterparts 
like T Tauri disks, the disk mass is orders of magnitude lower 
than the typical value found in T Tauri stars, and the scale height 
is independent on the mass of the central object. Moreover, we 
systematically compared the modeling results from Bayesian 
analysis between the ETBD and ETBD groups and found sim¬ 
ilar trends of flaring index as a function of SpT. The disk scale 
heights are comparable in both the high-mass and very low-mass 
BDs. Both groups feature a similar median disk mass and no 
clear trend is visible in the distribution, probably due to the un¬ 
certainty in translating the PACS photometry into disk mass, the 
detection bias, and the difference in evolutionary stage among 
the targets. Euture deeper far-IR surveys and ALMA observa¬ 
tions are required to improve constraints on the underlying mor¬ 
phology of the relationship between the stellar and disk proper¬ 
ties in the low-end stellar mass regime. 
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Appendix A: The Bayesian probabiiity distributions 
of seiected disk parameters for each modeied 
object 
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Fig. A.l. Bayesian probability distributions of p, //loo and Mdisk Fig. A.l. Continued, 
for each object. 
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